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Abstract — In this paper, a novel maximum Doppler spread 
estimation algorithm is presented for OFDM systems with the 
comb-type pilot pattern. The least squared estimated channel 
frequency responses (CFR's) on pilot tones are used to generate 
the auto-correlation matrices with/without a known lag, from 
which the time correlation function can be measured. The 
maximum Doppler spread is acquired by inverting the time 
correlation function. Since the noise term will bias the estimator, 
the estimated CFR's are projected onto the delay subspace of 
the channel to reduce the bias term as well as the computation 
complexity. Furthermore, the subspace tracking algorithm is 
adopted to automatically update the delay subspace. Simulation 
results demonstrate the proposed algorithm can quickly converge 
to the true values for a wide range of SNR's and Doppler spreads 
in Rayleigh fading channels. 

Index Terms — Doppler spread. Estimation, Subspace tracking, 
OFDM, Time correlation, Comb-type pilot. 

L Introduction 

In order to cope with various radio transmission scenarios, 
adaptive strategies, for example, adaptive modulation and 
coding and dynamic resource allocation, are widely employed 
by many orthogonal frequency division multiplexing (OFDM)- 
based wireless standai'ds, e.g., TAB, TVB, IEEE 802.11/16 
and 3GPP LTE [1]. Adaptive schemes automatically adjust the 
system configurations and transmission profiles according to 
some criteria to accommodate the varying radio environments. 

The maximum Doppler spread is one of the key parameters 
of criteria for adaptive strategies. It determines the fading 
rate of the radio channel and its reciprocal is a metric of the 
coherent time of the channel. With the knowledge of it, wire- 
less systems can change the depths of interleavers to reduce 
coding/decoding latencies, decrease unnecessary handoffs and 
adjust the rate of power control to reduce signalling overhead. 
More important, for OFDM systems, when the Doppler spread 
is comparable to the tone spacing, the orthogonality between 
tones would be corrupted, which would arise the inter-carrier 
interference (ICI) and consequently deteriorate the system 
performance. However, if the Doppler spread can be estimated. 
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it will facilitate the adaptive schedule/control algorithms to 
appropriately tune systems to alleviate ICI. 

Most of existing methods of estimating the maximum 
Doppler spread are categorized into two classes [2]: level 
crossing rate (LCR)-based and covariance (Cov)-based tech- 
niques. Since the algorithms reviewed in [2] were not specif- 
ically designed for OFDM systems, they did not exploit 
the special signal structure of OFDM systems. For OFDM 
systems, most literatures are Cov-based. [3] determined the 
maximum Doppler spread through estimating the smallest 
positive zero crossing point. In [4], Cai proposed to obtain 
the time auto-correlation function (TACF) in time domain by 
exploiting the cyclic prefix (CP) and its counterpart. However, 
Yucek [5] pointed out that for scalable OFDM systems whose 
CP sizes were varying over time, [4] was difficult to offer a 
sufficient estimation of TACF, which would degrade its estima- 
tion accuracy significantly. On the contrary, Yucek proposed 
to estimate the Doppler spread through the channel impulse 
responses (CIR's) which were estimated from the periodically 
inserted training symbols. 

Although the method in [5] seems to be more robust than in 
[4], its shortcomings are also evident. For example, in order to 
reduce system overheads, training symbols are arranged to be 
far from each other, and typically transmitted as preambles to 
facilitate the frame timing and carrier frequency synchroniza- 
tion. Once the duration of frame is longer than the coherent 
time of the channel, the maximum Doppler spread cannot be 
attained because TACF turns to irreversible. Moreover, for 
sparse training symbols, the converging speed of [5] would 
be very slow, which hinders its employment. 

In this paper, we propose to estimate TACF by exploiting 
the comb-type pilot tones [6] which are widely adopted in 
wireless standards. In order to reduce noise perturbation, the 
estimated channel frequency responses (CFR's) are projected 
onto the delay subspace [7] to obtain CIR's, and the subspace 
tracking algorithm [8] is adopted as well to track the drifting 
delay subspace. 

This paper is organized as follows. In Section the OFDM 
system and channel model are introduced. Then, the maximum 
Doppler spread estimation algorithm is presented in Section 
Hm Simulation results and analyses are provided in Section 
HVl Finally, Section IV] concludes the paper. 



II. System Model 

Consider an OFDM system with a bandwidth of BW ~ 
1/T Hz {T is the sampling period). N denotes the total 
number of tones, and a CP of length Lcp is inserted before 
each symbol to eliminate inter-block interference. Thus the 
whole symbol duration is Tg ^ (N + Lcp)T. In each OFDM 
symbol, P (< N) tones are used as pilots to assist channel 
estimation. In addition, optimal pilot pattern, i.e., equipowered 
and equispaced [9], is assumed. Pilot indexes are collected in 
the set Xp, i.e., Jp = {4> + p 9}, {p = 0, ...,P ~ 1), where 
(f> and 9 are the offset and interval, respectively. 

The discrete complex baseband representation of a multi- 
path CIR of length L can be described by [10] 

L-l 

h{n,T) = ^ ji{n)S [t -Ti) 

1=0 

where t; is the delay of the ^-th path, normalized by the 
sampling period T, and Jiin) is the corresponding complex 
amplitude. Due to the motion of users, 7/(n)'s are wide-sense 
stationary (WSS) narrowband complex Gaussian processes, 
and uncorrected with each other based on the assumption 
of uncorrelated scattering (US). In the sequel, P > L and 
P X 9 = N are assumed for determinability and simplicity, 
respectively. 

Furthermore, we assume the uniform scattering environment 
introduced by Clarke [11], thus 7i(n)'s have the identical 
normalized TACF jQ{2Trfdt) for all Vs, where fd is the 
maximum Doppler spread and Jo( ) is the zeroth order Bessel 
function of the first kind. Hence, the discrete TACF is 

rt4m)^afJoi27r\m\fdT) (1) 

where af is the power of the Z-th path. Additionally we assume 
the power of channel is normalized, i.e., ^f'^Q af ^ I. 

Assuming a sufficient CP, i.e., Lcp > L, the signal model 
in the time domain can be expressed as 

L-l 

ym{n) = ^ h,n{n, Ti)Xra{n - Ti) + (n) 

;=o 

where Xm{n) and ym{n) are the 7i-th samples of the m-th 
transmitted and received OFDM symbols, respectively, Wm {n) 
is the sample of additive white Gaussian noise (AWGN), 
i.e., E[wra{n)w':^{n + q)] = crld{q), and hm{n,Ti) is the 
corresponding sample of the time-varying CIR. Through some 
simple manipulations, the signal model in the frequency do- 
main is written as 

^ in ^^ni^m ^" ^^'m (2) 

where X,„,Y,„,W„i e qNxi jj^g jri-ih transmitted 
and received signal and noise vectors in the frequency do- 
main, respectively, and H„j e qNxn jj^g channel trans- 
fer matrix whose {y + /c, fc)-th element, i.e., \H.m\,j^k.k^ is 
TfY.n=oT.]:=o hm{n,n)e~^^<-^+^-^y^ , where k denotes 
subcarrier while v denotes Doppler spread. Apparently, as 
is non-diagonal, ICI is present. However, when the normalized 



maximum Doppler spread, i.e., fdTs, is less than 0.1, the 
signal-to-interference ratio (SIR) is over 17.8 dB [12], which 
enables us to discard non-diagonal elements of Hm with a 
negligible performance penalty. 

As the comb-type pilot pattern is adopted, only pilot tones, 
denoted as Yp;„i G C^^^, are extracted from Y„i. By 
approximating H„j to be diagonal, (|2]i is modified to 

Yp;„l ~ Xp;„lHp;„l + Wp;,„ (3) 

where Xp^m G C^^^ is a pre-known diagonal ma- 
trix, and Hp.,„ G C^^^ consists of diagonal elements 
of H„i. Hence, by denoting the instantaneous CFR as 
Hrain, k) = /^„,(n,T0e-^"2-'=-'/^, we have [Hp.„]p = 

7f'^n=o ^rni'iT'T'p + P X 9). Denote the Fourier transform 
matrix on the pilot tones as Fp G C^^^, that is, [Fp]p.„ = 
i^-j2n{4'+Pxe)n/N^ jjjgjj Wp..n^ = Fpw,„, whcrc w™ = 

[w,n{0), . . .,Wm{N - SO, E[Wp,„,W^.J = allp. 

III. Maximum Doppler Spread Estimation 

At the receiver, the least-squared (LS) channel estimation 
on pilot tones is carried out firstly, i.e., 

Hp;m = Xp.-'j^jYp;,,! = Hp;„l + Vp:„J (4) 

where flp-m G C^^^ is the estimated CFR, and Vp-m G 
C^^^ is the noise vector expressed as \p:m — Xp.^^Wp:m, 
hence, Vp;,„ ~ CA/'(0, cr^Ip) when Xp.„jXp;,„ — Ip for 
PSK modulated pilot tones with equal power 

In the following, we will introduce a method of estimating 
the maximum Doppler spread based on TACF measured from 
significant paths of the channel obtained by projecting the LS 
estimated CFR onto the delay subspace. 

A. Measurement of the Time Auto-Correlation Function 

First, by defining the Fourier transform matrix as Fp ,- G 
CP^^ with [Fp,^]p,, = g-j2,r(0+pxe)r,/jv^ Yip.r,,. Can be 
expressed as 

Hp:m = ^ X! Hp;„i(n) = X! Fp,Thm("-) 

ri=0 ri=0 

where Hp;m and hm(n) are CFR and instantaneous CIR 
vectors, respectively. Regardless of noise, the 0-lag auto- 
correlation matrix of Hp:„, is 

R-Hp(O) = E [Hp;mHp.„] 

= ^Y.Y. Fp,rA™(n, q)^% (5) 

n=0 g=0 

where A„i{n,q) = E[hm{n,)h^{q)], and based on the as- 
sumption of WSSUS and Clarke model, its {li, Z2)-th element 
is [A„i{n,q)]i^ i^ = crfjt{n - q)S{li - h), where rt{n) is 
the normalized TACF, hence Am{n,q) is diagonal. Denoting 
D — diag{af), Z = 0, . . . , L — 1, we have 

Ara{n,q) ^ rt{n - q)T) (6) 



Substitute ^ into ^ and with some manipulations 

Rhp(O) = ^(0)Fp,,DF^^, (7) 

Af-lAf-l 

^(o) = ]^EE (8) 

n=Q q=0 

Similarly, the /?-lag auto-correlation matrix of Hp-m iP > 0), 
defined as Rhp(/3) — -E[Hp:.,„+/3Hp.„J, can be written as 

Rhp(/3) = e(/3)Fp,.DF«, (9) 

JV-lJV-l 

^(^) = ]^ E E - 9 + + ^-p)^) ^10) 

Then, with (|7]i and (|9]l, we have 



|Rhp(/?)||| 

2 
F 



(11) 



IRhp(0)||^p e(0) 

where 1 1 • 1 1 f denotes the Frobenius norm. When the normalized 
Doppler spread JdTs < 0.1, referring to ([TJ, we can make an 
approximation (which we will examine later) as 

r, « Jo{27Tp{N + L,p)fdT) (12) 

Since when /3{N+Lcp)fdT = /S/dT^ < 0.38, Jo(-) is positive 
and reversible, meanwhile, in order to hold the orthogonality 
between subcarriers, fdTs is usually smaller than 0.1, thus 
/3 < 3 is the feasible range. Then fd can be estimated by 



fd — 



(13) 



Now we consider the effect of noise. When noise is present, 
O and (|9]l are rewritten into 

Rhp (0) = e(0)Fp,,DFf , + allp (14) 

Rhp(/3) = ^(/3)Fp,,DFf_, (15) 

Correspondingly, (fTTT l changes into 



77 : 



|Rhp(/3)II| 



IRhj 



where p is defined as 

P = 



I ||Fp,.DF|^^J||, 
As pilot tones are equispaced, Fp^Fp^ = PIl, then 



(16) 



(17) 



|Fp,DFf,J|| = P^Va 



L-l 
'2 \ " „4 



1=0 



therefore, ( fTTl ) is 



P = 



(18) 



Note Ef=o' '^f < (EtV '^f we have 

1 

y/PxSNR 

where SNR = cr~^ for the normalized power of the channel 
and pilot tones. 



(19) 



B. Improving the Estimation Accuracy by the Delay Space 

Although the maximum Doppler spread can be evaluated 
from (fTST i and ( fTsT l, the effect of noise will bias the result 
of estimation heavily when P is small and SNR is low. 
On the other hand, when P is large, the effect of noise is 
negligible, but the sizes of Rhp (0) and Rhp (/?) turn to be 
so large that evaluating their Frobenius norms would require a 
large amount of calculations, which hinders the employment 
of this method in real applications. Therefore, we introduce 
the delay subspace onto which the auto-correlation matrices 
are projected to reduce the effect of noise as well as the 
computation complexity. 

Firstly, the eigenvalue decomposition (EVD) is performed 

Rhp (0) = [e(0)A + allp] Uf (20) 
Rhp(/3) = U,K(/3)A]Uf (21) 

Since the number of channel taps is L, rank{A) ^ L < P, 
the last P ~ L eigenvalues of Rhp (0) and Rhp (/?) are cr^ 
and 0, respectively. Once L is available, Ur can be divided 
into two parts named as the "signal" and "noise" subspaces, 
respectively, i.e., U,- = [Ui-,s, U-r,n], where IJr.s G C^^^, 
and so does A, i.e., A = diag{As, Op-l), where A^ G C^^^. 
Hence, 



-r,s-l-Hp(0)Ur 



e(0)A, + allL 

,Hp(/3)u,,. = m^s 

Based on ( |22] | and ( l23T l. ( fTSI l can be refined as 



Uf.R: 



V = 



||U^,Rhp(/3)U,. 



sIlF 



\ ||U^,Rhp(0)U.^ 
where pr is defined as 



Pr 



Pr = 




(22) 

(23) 



(24) 



(25) 



From ([IliOdaOlldlB, 



|Fp.,DF 



H 

P.T 



= 1 1 A., 



1^ 



Hence, comparing (l25l l with ( [TtI i. the bias term is reduced 
because 

Pr I L 



< 1 



(26) 



Actually, in the real circumstance, the number of significant 
taps of wireless channels is far less than of pilot tones, thereby 
projecting auto-correlation matrices onto the delay subspace, 
like (|22] | and ( |23] |, can effectively reduce the bias term and 
ease the calculation of 77. 

C. Tracking the Delay Subspace - the Proposed Algorithm 

When the user is moving, the tap delays of the channel, i.e., 
T;'s, are slowly drifting [13] [7], which causes Fp;t to vary 
and so does XJt.s - To accommodate this variation, the subspace 
tracking algorithm is adopted to automatically update the delay 
subspace. In addition, if the number of significant taps of 



the channel is unknown, minimum description length (MDL) 
[14] is employed to estimate it. The proposed algorithm is 
summarized as follows. 

Initialize: (n = 0) 

Qo(0) = Q^(0) = [lL,„,OL.p_i„]^ 
Ao(0) = A^(0) = Op,i„ 
Co(0) =0^(0) =Il„. 
Run : ( n = 71 + 1 ) 
Input: Hp(n) 

1) Updating for the 0-lag auto-correlation 
matrix : 

Zo(n) = Qo(n- l)Hp(n) 

Ao(n) = QAo(n - l)Co(n - 1) + (1 - a)Hp(n)Zo(n)" 
Ao(n) = Qo(n)Ro(n) : QR-factorizat±on 



TABLE I 
A Table of Values of ( fT2t 



Co(n) = Qo(r 



'Qo(n) 



L(n) = MDL (dmg(Ro(n))) 

2) Updating for the /3-lag auto-correlation 
matrix : 

Z0{n) =Q^(n-l)Hp(n-/3) 

A^(n) = aA0{n - l)C^(n - 1) + (1 - Q)Hp(n)Z/3(n)^ 
A/3(n) = Q/3(n)R;3(n) : QR-factorizat±on 

Ci3{n) =Qp{n-l)"Q0{n) 

3) Estimating 77 according to (1241 ) : 



^ Efir'l[Ro(n)],J 



4) Estimating /d according to (1131) . 

Remark: a is a positive exponential 
forgetting factor close to 1 . Lm is the 
maximum rank to be tested. MDL{-) denotes 
the MDL detector and diag(R.) denotes the 
diagonal elements of R. In the simulation, 
we set Q = 0.995 and Lm ~ W . 



D. Other Considerations 

In this subsection, several further discussions about the 
proposed algorithm are presented. 

First, numerical results are shown in Table U to examine 
( fnt when N = 512, L^p = 64 and /3 = 1. From Table H] we 
find ( fT2] i is a good approximation when PfdTg is small. It is 

worth noting that Jo{2t: (3 fdTs) < 77, hence fd > ^^r^ 



m 



other words, (fTZt tends to over-estimate the maximum Doppler 
spread a bit. 

Then we compare the computation complexity of the pro- 
posed algorithm with (fTSI l. The computation complexity of the 
subspace tracking is 0{P x L^) [8]. Since it takes a dominant 
proportion of the total number of instructions required by 
the proposed algorithm, we use it instead. Meanwhile, the 
computation complexity of (fT6] l. which directly computes the 
Frobenius norm of P x P matrices, is 0{P^). Apparently, 
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Fig. 1. Performance comparison for the CP-based [4], Frobenius-norm- 
based )16t and subspace-based j24t algorithms when a 20ms frame is used 
and /3 = 1. 



when P > L^, which is usually the case for sparse multipath 
channels, the proposed algorithm can reduce the computation 
complexity considerably. 

IV. Simulation Results 

The performance of the proposed algorithm is evaluated for 
an OFDM system with BW = 5 MHz (T = l/BW = 200 
ns), iV = 512, Lcp = 64 and P = 64. ITU Vehicular A Chan- 
nels [15] is adopted, which consists of six individually faded 
taps with relative delays as [0, 310, 710, 1090, 1730, 2510] ns 
and average power as [0, —1, —9, —10, —15, —20] dB. The 
classic Doppler spectrum, i.e., Jakes' spectrum [10], is applied 
to generate the Rayleigh fading channel. 

In FiglU a CP-based algorithm reported in [4] and ( [T6b . 
which is based on the Frobenius norm, are compared with the 
proposed subspace-based algorithm (l24l l with P ~ 1. A 20ms 
frame including 194 OFDM symbols is used to generate the 
statistics, fd = 200, 400 and 600 Hz are tested under a range 
of SNR's, respectively. Apparently, the CP-based algorithm 
fails for all /^'s when the SNR is below 20 dB, meanwhile 
( fT6] l and (OH work very well for all fd's and SNR's but with 
a moderate positive bias for SNR = 5 dB. In fact, when 
SNR = 5 dB, resorting to iT% and (l26b . the lower bound 
of the bias terms p and pr are 0.0395 for (O and 0.0099 
for (l24l) . respectively. And when fd = 600 Hz, according to 
IHl and (Uni), ^(0) = 0.9938 and f (/3) = 0.9476. Thus, the 
relative errors of r] are 0.0007 and 0.0000 for (US)) and (|24|, 
respectively, which almost have no effect on the estimation of 
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Fig. 2. Performance comparison between the Frobenius-norm-based )16t and 
subspace-based I24t algorithms when = 400 Hz and /3 = 1. 



Fig. 3. The convergence of the proposed subspace-tracking-based algorithm 
when SNR = 15 dB and /3 = 1. 



the maximum Doppler spread. Therefore, ( fTSI l and (|24] | show 
almost the same performance when SNR's are above 5 dB. 

Fig|2] shows the performance comparison between the 
Frobenus-norm-based ( fT6] l and subspace-based ( |24] | in the 
low SNR regime, specifically, below 5 dB, to emphasize the 
capability of noise depression of the latter. Different frame 
durations are used to obtain TACF. From the figure we can find 
(l24l l outperforms (fTSI l for all the SNR's and frame durations, 
although both of them over-estimate the maximum Doppler 
spread due to the non-negligible noise bias term in the low 
SNR regime, which is also the reason why increasing the 
length of observation record does not help to decrease the 
bias in this regime. 

The convergence of the proposed subspace-tracking-based 
algorithm is shown in Fig|3] Three different maximum Doppler 
spreads are tested, i.e., f^i — 200, 400 and 600 Hz, when 
SNR = 15 dB. The curves are plotted from the thirtieth 
OFDM symbol. It is observed from the plot that all the 
three curves converge to their true values after numbers of 
OFDM symbols and, further, the higher the Doppler spread 
is, the faster the curve converges. This is due to the subspace 
is updating faster when the Doppler spread is higher. In 
addition, after converging, the estimated maximum Doppler 
spread is fluctuating around its true value in a small range, 
hence additional time-averaging can be employed to smooth 
the curve. 

V. Conclusions 

In this paper, we propose a subspace-tracking-based maxi- 
mum Doppler spread estimation algorithm which is applicable 
to OFDM systems with the comb-type pilot pattern. It enjoys 
three main advantages: i) alleviating the noise term; ii) re- 
ducing the computation complexity; iii) tracking the drifting 
delay subspace. Through simulations, the performance of the 
proposed algorithm is demonstrated to outperform the CP- 
based algorithm [4]. Moreover, since the proposed algorithm 
is based on the subspace tracking, it can be easily integrated 



into the channel estimator equipped with the subspace tracker 
[7] [16], which lends a broad application promise to it. 
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